;+
; NAME:
;       EXTINCTION_COLUMN_MAPS
;
; PURPOSE:
;       Calculate N(H2) based on the 8-um GLIMPSE maps and the
;       IRDC-Modeling Galactic 8-um emission model.
;
; CATEGORY:
;       IRDC modeling routine
;
; CALLING SEQUENCE:
;       EXTINCTION_COLUMN_MAPS, struct, dist [,VLSR=vlsr][,/NEAR][,/FAR]
;
; INPUTS:
;       STRUCT -- The 1-element BGPS catalog structure for a single
;                 BGPS object (see READ_BGPS_CSV.pro).
;       DIST   -- Distace to source [kpc].  Optionally, may supply
;                 VLSR to calculate kinematic distance.
;
; OPTIONAL INPUTS:
;       VLSR   -- The VLSR of the object for the calculation of
;                 kinematic distance.  Will use NEAR kinematic
;                 distance, unless /FAR keyword specified.
;
; KEYWORD PARAMETERS:
;       NEAR   -- Use the NEAR kinematic distance [default].
;       FAR    -- Use the FAR kinematic distance.
;       RV31   -- Use the Bohlin-Savage conversion for R_V = 3.1 dust
;                 (duffuse ISM) [Default: use conversion for R_V = [5.5]
;       PLOT   -- Create an .EPS plot in ./local/output/
;
; OUTPUTS:
;       NONE    [Writes maps to .FITS files.]
;
; OPTIONAL OUTPUTS:
;       Uknown
;
; MODIFICATION HISTORY:
;
;       Created:  08/20/12, TPEB -- Initial Version.
;
;-

PRO EXTINCTION_COLUMN_MAPS, s, dist, VLSR=vlsr, NEAR=near, FAR=far, RV31=rv31, $
                            PLOT=plot
  
  IF n_elements(s) EQ 0 THEN message,'Requires a Bolocat structure element.'
  plot = KEYWORD_SET(plot)
  
  scnum = string(s.cnum,format="(I4)")
  
  ;; Check that the UBC_GLIMPSE_PROC postage stamps exist
  DIR = './local/ubc_glimpse_proc/'
  IF ~ FILE_TEST(DIR+'bgps'+scnum+'.fits') THEN BEGIN
     message,'BGPS #'+scnum+' is not covered by GLIMPSE -- '+$
             'no column file written.',/cont
     RETURN
  ENDIF
  
  ;; Check the distance / velocity
  IF n_elements(dist) + n_elements(vlsr) EQ 0 THEN BEGIN
     message,'Routine requires a either a user-supplied distance estimate '+$
             'or a VLSR to calculate a kinematic distance!',/cont
     RETURN
  ENDIF
  
  ;; Get distance in pc
  IF n_elements(dist) EQ 0 THEN BEGIN
     dist = (KEYWORD_SET(near) || ~KEYWORD_SET(far)) ? $
            KDIST(s.glon_peak,s.glat_peak,vlsr,/NEAR) : $
            KDIST(s.glon_peak,s.glat_peak,vlsr,/FAR)
  ENDIF ELSE BEGIN
     ;; Check that supplied distance is in kpc (Guess < 100 in kpc).
     IF dist LE 100. THEN BEGIN
        message,'Converting distance in kpc to pc...',/inf
        dist *= 1.d3
     ENDIF
  ENDELSE
  message,'Distance = '+string(dist,format="(F0.1)")+' pc',/inf
  
  ;; Read in postage stamp images
  imir = readfits(DIR+'Imir'+scnum+'.fits',imir_hdr,/silent)
  irac = readfits(DIR+'smooth30arc'+scnum+'.fits',irac_hdr,/silent)
  IF plot THEN smoo = readfits(DIR+'bsmoo'+scnum+'.fits',bsmo_hdr,/silent)
  
  ;; Calculate C, f_fore, & N(H2)
  c = (1. - irac / imir) > 0.
  ffore = ffore_lbd(s.glon_peak,s.glat_peak,dist,/NONNATIVE)
  
  print,dist
  
  column0 = KEYWORD_SET(rv31) ? 1.915d22 : 2.816d22
  nh2 = (column0 * (-alog(1.d - c / (1.d - ffore))))>0.
  
  ;; Prep the column map for writing to FITS file
  sxaddhist, '==============================', irac_hdr, /COMMENT
  hist = ['EXTINCTION_COLUMN_MAPS: '+systime(0),$
          'Created 8-um Extinction Column Density Map']
  sxaddhist, hist, irac_hdr
  sxaddpar, irac_hdr, 'OBJNAME',s.name,' BGPS Source Name'
  sxaddpar, irac_hdr, 'BGPSCNUM',s.cnum,' BGPS Catalog Number'
  sxaddpar, irac_hdr, 'SMOOIRAC','smooth30arc'+scnum+'.fits',$
            ' Smoothed (30") GLIMPSE cutout'
  sxaddpar, irac_hdr, 'IMIR_EST','Imir'+scnum+'.fits',$
            ' Estimate of total mid-infrared, I_MIR'
  sxaddpar, irac_hdr, 'IMIR_AVE',mean(imir),' [MJy/sr] Mean I_MIR for tile',$
            FORMAT="F0.1"
  wfn = 'ir_column'+scnum+'.fits'
  sxaddpar,irac_hdr,'FILENAME',wfn
  sxaddpar,irac_hdr,'BUNIT','cm^-2'
  sxdelpar,irac_hdr,['JY2DN','ZODY_AVE','GAIN','FLUXCONV']
  rvused = KEYWORD_SET(rv31) ? 3.1 : 5.5
  sxaddpar,irac_hdr,'RV_USED',rvused,' Selective extinction ratio (R_V)',$
           FORMAT="F0.1"
  sxaddpar,irac_hdr,'COLUMN_0',column0,' [cm^-2] Column prefactor (Eqn. 13)',$
           FORMAT="E0.3"
  sxaddpar,irac_hdr,'DISTANCE',dist/1.d3,' [kpc] Heliocentric distance',$
           FORMAT="F0.2"
  sxaddpar,irac_hdr,'F_FORE',ffore,' Foreground fraction from model',$
           FORMAT="F0.3"
  
  writefits,DIR+wfn,nh2,irac_hdr
  
  
  
  
  ;;===========================================================
  ;; Plotting junk....
  IF plot THEN BEGIN
     myps,'./local/output/ir_column'+scnum+'.eps',xsize=11,ysize=10
     multiplot,[2,2],gap=0.04,mtitle='BGPS #'+scnum
     
     ;;=========
     ;; Panel #1 -- Smoothed IRAC Band 4 + BGPS contours
     undefine,pr
     plotmap,irac,irac_hdr,ct=0,range=pr,charsize=1.0,min_dpi=50
     cgColorbar,/vert,/right,position=[0.94,0.56,0.96,0.93],range=pr,$
                title='GLIMPSE Intensity [MJy sr!u-1!n]',charsize=1.0
     plot_wcs_axes,bsmo_hdr,10,lonarr=bxc,latarr=byc
     cgLoadct,40,/silent
     contour,smoo,bxc,byc,levels=[.025,.05,0.1,0.25,0.5],$
             c_colors=[1,2,3,4,5]*50,/overplot,thick=2
     multiplot
     
     ;;=========
     ;; Panel #2 -- I_MIR estimate
     plotmap,imir,imir_hdr,ct=0,range=pr,charsize=1.0;,min_dpi=50
     multiplot
     
     ;;=========
     ;; Panel #3 -- Contrast calculated from above
     plotmap,c,imir_hdr,ct=13,range=[0,1],charsize=1.0,min_dpi=50
     cgColorbar,position=[0.17,0.04,0.49,0.06],range=[0,1],$
                title='Contrast',charsize=0.7
     multiplot
     
     ;;=========
     ;; Panel #4 -- Column Density, calculated from #3 + ffore
     cpr = [0,3.d22]
     plotmap,nh2,imir_hdr,ct=41,range=cpr,charsize=1.,min_dpi=50
     cgColorbar,position=[0.59,0.04,0.92,0.06],range=cpr,$
                title='N(H!d2!n) [cm!u-2!n]',charsize=0.7
     
     myps,/done,/mp
  ENDIF
  ;; END Plotting junk....
  ;;===========================================================
  
  
END
